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SUMMARY 

A unique modification to the NASTRAN solution sequence for transient 
analysis with direct time integration (COSMIC NASTRAN rigid format 9) has been 
developed and incorporated into a DMAP alter. This DMAP alter calculates the 
buckling stability of a dynamically loaded structure, and is used to predict 
the onset of structural buckling under stress-wave loading conditions. The 
modified solution sequence incorporates the linear buckling analysis capability 
(rigid format 5) of NASTRAN into the existing Transient solution rigid format 
in such a way as to provide a time dependent eigensolution which is used to 
assess the buckling stability of the structure as it responds to the impulsive 
load. As a demonstration of the validity of this modified solution procedure, 
the dynamic buckling of a prismatic bar subjected to an impulsive longitudinal 
compression is analyzed and compared to the known theoretical solution. In 
addition, a dynamic buckling analysis is performed for the analytically less 
tractable problem of the localized dynamic buckling of an initially flawed com- 
posite laminate under transverse impact loading. The addition of this DMAP 
alter to the transient solution sequence in NASTRAN facilitates the computa- 
tional prediction of both the time at which the onset of dynamic buckling 
occurs in an impulsively loaded structure, and the dynamic buckling mode 
shapes of that structure. 


INTRODUCTION 

Composite laminates that are subjected to static, dynamic, or fatigue 
loading are known to undergo delamination, or debonding, between the laminated 
plies of which they are composed. Delamination causes a significant loss 
stiffness and strength, and can considerably reduce the structural integrity 
of a laminate. Once this damage has occurred, a compressive stress near the 
delamination can induce local buckling of the delaminated plies. This buck- 
ling may then cause further extension of the delamination and progressive wea- 
kening of the laminate. In lieu of actual experimental testing, the ability 
to computationally predict the onset of delamination buckling is necessary for 
evaluating the durability of many composite structures. 

The delamination buckling phenomenon has been observed experimentally 
under both static and fatigue loading conditions (Refs. 1 to 4), and several 
analytical and numerical methods have been proposed (Refs. 5 and 6) to model 
this damage mechanism. Finite-element approaches (Refs. 7 to 9) have been 
used as the basis for these analyses, but no comparable numerical methods exist 
to analyze delimination buckling which occurs as a result of an impulsively 
applied load. That is the topic of this paper. 
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Exper i-mentaf- observations of dynamic delamination buckling in transversely 
impacted laminates were reported earlier (Refs. 10 to 12), using high-speed 
photography and simultaneous strain measurements of transversely impacted lami- 
nates. A related numerical analysis (Ref. 10) indicated that the buckling 
behavior must be accounted for in the computational model in order to accu- 
rately assess the damage tolerance capability of the laminate. This motivated 
the present development of a NASTRAN DMAP alter analysis procedure that can 
be used to computationally predict the onset of buckling instability under 
transient stress-wave loading. 

The objectives of this paper are, therefore, (1) to outline the dynamic 
buckling analysis computational procedure and its implementation into the DMAP 
alter sequence (2) demonstrate the validity of the dynamic buckling analysis 
procedure by analyzing a simple one-dimensional example problem with a known 
solution, and (3) apply the dynamic buckling analysis to the analytically less 
tractable problem of the localized dynamic buckling of an initially flawed 
composite laminate under transverse impact loading. 

The NASTRAN transient solution sequence, when modified as indicated in 
the following section, provides a new computational tool that can be used to 
predict both the time at which the onset of dynamic buckling occurs and the 
dynamic buckling mode shapes of an impulsively loaded structure. 


Dynamic Buckling Analysis 

Linear buckling analysis requires solution of the eigenvalues problem: 

[[K] + XCK a ] {<*>} - o] (1) 


where 

[K] structural stiffness matrix; 

[K a ] stress stiffness matrix 

X, {<t>} denote the associated eigenvalue and eigenvector 

In terms of the buckling analysis, the eigenvector {<j>} represents the 
buckling mode shape, and the associated eigenvalue X indicates the multiple 
of [K a ] needed to make equation (1) singular, that is, to cause buckling. In 
a one-dimensional column buckling problem, each scalar eigenvalue satisfying 
equation (1) physically represents the nondimensional ratio: 



where a is the compressive stress in the column, A is the cross-sectional 
area, and P* is the buckling load. If the eigenvalue has the critical value 
of unity (oA = P*), buckling in the associated mode occurs. 

In the dynamic case, the terms of CK a ] in Eq. (1) vary with time as the 
stress waves propagate through the structure. The eigensolution of (1) then 
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becomes time dependent, and can be used to track the buckling stability as a 
function of time. Figure 1 is a simplified representation of a modified 
direct-time integration solution sequence in which the updated stress stiffness 
matrix is formed after each time step At, and the associated eigenvalue prob- 
lem in equation (1) is solved. The eigenvalue is now a function of time, and 
it indicates the onset of buckling when it reaches the critical value of unity. 
Figure 2 is the DMA P alter which incorporates this dynamic buckling algorithm 
into the existing transient solution sequence. 


DMAP Procedure 

The functions of the DMAP statements shown in Fig. 2 are summarized 
here. In line 2 the number of columns in the UPV matrix is determined. This 
matrix contains the displacement, velocity and acceleration vectors for each 
degree of freedom at each time step. Lines 2 through 16 follow the Bubble 
Algorithm approach of Ref. 13. The DMI column matrices TIPI and BAS1 from the 
Bulk Data deck, each initially sized to contain more rows than columns in the 
UPV matrix, are used to form two new column matrixes, MNTRJ and BOOTI. The 
number of rows in each of these matrices is equal to the number of columns in 
the UPV matrix. The monitor matrix MNTRJ initially contains unity in the 
first row and zero in the remaining rows. The BOOTI matrix always contains 
unity in the last row and zero in the remaining rows. 

Having determined the size of the partitioning matrices, the eigenvalue 
extraction data is determined in line 19 and the buckling calculations are now 
performed. At the beginning of each pass through the RAALOOP, corresponding 
to each integration time step of the requested output, the current column 
position is compared with the number of columns in the UPV matrix, lines 25 
through 27, ending the loop at the end of the available data. Continuing 
within the loop the unity value of the MNTRJ matrix is advanced three rows, 
lines 28 through 31, pointing to the location of the current displacement 
vector in the UPV matrix. The MNTRJ matrix is used to partition the UPV 
matrix, line 32, stripping the column containing the displacements. These 
displacements are used in the DSMG1 module, line 33, to form the time-varying 
global differential stiffness matrix, KDGG. The reduced differential stiff- 
ness matrix, KDAA, is then formed by eliminating the restrained and dependent 
degrees of freedom, line 35 through 45, and in line 47 this matrix is multi- 
plied by negative one, forming the KDAAM matrix. The stiffness matrices KAA 
and KDAAM are then used in the READ module, line 48, to solve for the eigen- 
values and eigenvectors for each integration time step initially requested for 
output. 

The eigenvalue for each time step is printed by line 52. Optionally, 
lines 53 and 54 may be used to print eigenvalues and eigenvalue extraction 
data. Line 58 may be used to print eigenvectors. The RAALOOP is ended at 
line 64. 

The computationally intensive nature of this analysis can be made more 
efficient by slightly modifying the DAMP procedure. A promising method is to 
perform the buckling analysis at specified time intervals in the transient 
solution sequence rather than after every time step, as is done here. The 
length of the time interval can be progressively decreased as the eigenvalue 
begins to change more rapidly, or as the critical value of unity is approached. 
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This technique will significantly reduce the number of individual buckling 
analyses performed, and hence will result in a more computationally efficient 
algorithm. 


Example Problem 

In order to establish the validity of this analysis procedure, a simple 
problem with a known solution, as given in Ref. 14, was analyzed. The propaga- 
tion of a longitudinal compressive pulse in a long prismatic bar, shown in 
Fig. 3, was modelled. 

Assuming a one-inch diameter aluminum bar of uniform circular cross sec- 
tion the elastic and geometric constants are: 


E = 10 x 10 6 psi 

(3) 

4 4 

I * irr = ir in. 

4 64 

(4) 

2 2 
A = irr = ir in. 

4 

(5) 

„ r.,in-4 lb • s 2 
p = 2.5x10 4 

in. 

(6) 

L = 100 in. 

(6) 


where E is the Young's Modulus, I is the area moment of inertia, A is the 
cross-sectional area, p is the mass density, and L is the length of the bar. 

The lowest buckling load is given by (Ref. 15): 

P* = thl = 121 lb (7) 


As shown in Fig. 3, the applied load is identical to the static buckling load 
in Eq. (7) . 

Using the above material constants, the bar wave velocity is given by 
(Ref. 14): 


C 0 .yr . 200.000 if, (8) 

so the time for the longitudinal compression wave to travel from the impact 
point to the distal end of the bar is 


t = J=- = 500 ps 
0 C o 


(9) 
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A NASTRAN model consisting of ten rod elements, for a total of ten uncon- 
strained axial degrees of freedom, was used to model the longitudinal impact 
of the bar. The integration time step was taken as 

At = 4 TOC" = 12-5 pS (10) 

o 

to insure a numerically converged solution. The propagation of the compres- 
sion wave from the point of impact to the clamped end of the bar is depicted 
in Figs . 4(a) and (b) . 

The compressive pulse, traveling at a speed C 0 , reaches the complete 
length of the bar at time t 0 (500 ys). Because the distal end of the bar is 
held fixed, the incident compressive pulse reflects (Ref. 15) as a pulse of 
the same sign (compressi ve) which superimposes on the existing uniform compres- 
sive stress in the bar. Figures 4(c) and (d) depict the progression of the 
reflected pulse, traveling at a speed C 0 , back to the proximal end of bar, 
effectively doubling the compressive load supported by the bar. Reflecting 
from the proximal (free) end as a pulse of opposite sign (tensile) which super- 
imposes on the existing compressive stress, the bar returns to its original 
fully stressed state at time 3t 0 , (1500 ys) as shown in Figs. 4(e) and (f). 
Finally, in Figs. 4(g) and (h), the tensile pulse reflects as a tensile pulse 
from the fixed end which temporarily cancels the uniform compression at time 
4t 0 (2000 ys), leaving the bar instantaneously unstressed. The stress states 
depicted in Figs. 4(1) and (j), for all practical purposes identical to those 
in Figs. 4(a) and (b), indicate that, assuming no damping exists, the above 
cycle will repeat itself indefinitely. 

The corresponding time dependence of the lowest eigenvalue is shown in 
Fig. 5. The critical value of 1.0 is reached at times t 0 , 3t 0 , 5t 0 , 7t 0 ,. . . 
(500, 1500, 2500, 3500 ys,. . .); and whenever the bar supports a uniform com- 
pressive stress corresponding to its buckling load. Similarly, the eigenvalue 
reaches to its lower limit of 0.5 at times 2tg, 6t 0 , 1 0 t 0 , . . . (1000, 3000, 
5000 ys,. . .); and whenever the stress state is double that of the buckling 
load. The eigenvalue becomes large (theoretically infinite) at time 0, 4t 0 , 
8t 0 ,. . . (0, 2000, 4000, 6000 ys , . . .) ; and whenever the bar is unstressed. 

Superimposed on the finite element results in Fig. 4 is the theoretical 
1-D solution, assuming the stress wave propagates nondi speri vely at a constant 
speed C 0 and reflects from the boundaries as described above. Good agree- 
ment exists between the two solutions, even when relatively few finite ele- 
ments are used to model the bar. The time behavior of the lowest eigenvalue, 
shown in Fig. 5, can be interpreted directly in terms of the transient stress 
distribution in Fig. 4. Since the applied compressive load is exactly equal 
to the first static buckling load in Eq. (7), and no strain-rate dependence 
was assumed in the finite element model, buckling is predicted whenever the 
bar is uniformly stressed with its critical static buckling stress, which 
occurs at odd multiples of t 0 , as shown in Fig. 4. 

In a practical application, the above analysis is valid only until the 
onset of buckling occurs, since no post-buckling behavior has yet been included 
in the finite element model. The time itegration was extended in the example 
problem only to physically interpret the results of the dynamic buckling 
analysis . 
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Dynamic Delamination Buckling 


The example problem could have been solved without the use of a finite 
element analysis because of the simple non-di spersi ve nature of the longitudi- 
nal wave propagation. However, the propagation of flexural waves in beam-like 
structures is dispersive by nature, and as such would pose a formidable chal- 
lenge without the use of some type of computational simulation. In Ref. 11, 
experimental measurements of delamination duckling in graphite/epoxy composite 
laminates were reported. The beam-like experimental specimens had simulated 
delaminations (ply disbonds) embedded in them during the fabrication process. 
They were held clamped at both ends and impacted transversely, as depicted 
schematically in Fig. 6. The subsequent flexure-induced local buckling of the 
delamination was recorded using strain gages and high speed photography. A 
finite element model of the initially flawed experimental specimen is used 
here to verify that the dynamic delamination buckling phenomenon can be pre- 
dicted using computational simulation. Figure 6 shows the geometry and loading 
conditions for the initially flawed composite laminate subjected to a trans- 
verse impact. The finite element discretization of this laminate near the 
embedded flaw is shown schematically in Fig. 7. The layered structure of the 
composite laminate is represented by layers of shell elements. Multipoint con- 
straints are imposed on the degrees of freedom between neighboring nodal points 
in the thickness direction such that simple beam bending displacements are 
enforced; that is, plane sections remain plane and no strain exists in the 
thickness direction. These constraints are removed in the delaminated region 
to allow the delaminated plies to separate from the main laminate when a local 
compression occurs in that area, as shown in Fig. 7. More complete details of 
the finite element modeling procedure are given in Ref. 12. 

The progression of the flexural waves out from the central impact point 
to the boundaries of the laminate are shown in Fig. 8. As the disturbance 
passes through the flawed region at 100 to 150 ps after impact, the delaminated 
ligament separates from the laminate and begins to support a compressive longi- 
tudinal stress which increases in magnitude until it causes a local buckling 
of the delamination. The eigenvalue behavior and corresponding buckling mode 
are shown in Fig. 9. As the laminate deforms under the applied load, the 
eigenvalue decreases monotoni cal ly in magnitude until it reaches the critical 
value of unity, indicating the onset of buckling at approximately 190 ps from 
impact. The corresponding buckling mode shape is also depicted in the figure. 

These results correspond closely with experimental observations. Both 
the buckling mode shape and the time at which buckling occurs are in good 
agreement with measurements taken from high speed photographs. A detailed com- 
parison of finite element results and experimental measurements is given in 
Ref. 11. 


CONCLUSIONS 

A dynamic delamination buckling analysis procedure has been incorporated, 
in the form of a DMAP alter, into the transient analysis rigid format of 
NASTRAN. With this enhancement, NASTRAN can be used to calculate the time at 
which dynamic buckling occurs and the buckling mode shape of a structure sub- 
jected to dynamic loading. Comparison of the calculated results with a known 
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solution supports the validity of the analysis. Application of the dynamic 
buckling analysis to the more complex problem of transverse impact of beam- 
like laminate was demonstrated, and the results phenomenologically duplicated 
those reported in earlier experiments. 
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FIGURE 1. - DYNAMIC BUCKLING ANALYSIS SOLUTION SEQUENCE. 
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FIGURE 2. - DYNAMIC BUCKLING DMAP ALTER FOR TRANSIENT SOLUTION SEQUENCE. 
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FIGURE 5. - TIME DEPENDENCE OF FIRST EIGENVALUE IN 
EXAMPLE PROBLEM. 




FIGURE 7. - FINITE ELEMENT REPRESENTATION 
OF COMPOSITE LAMINATE SHOWING LOCAL BUCK- 
LING OF THE DELAMINATED PLIES. 
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FIGURE 9. - TIME DEPENDENCE OF LOWEST EIGENVALUE FOR COM- 
POSITE LAMINATE SUBJECTED TO TRANSVERSE IMPACT. 
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